ode — ordinary differential equation solver


\begin{rail}
ODE : 'ode' '(' Function ',' Matrix ',' Scalar ',' Scalar ( ',' (Scalar (',' ( Scalar (',' (Scalar ?)?)?)?)?)?) ')' ;
\end{rail}
ode is a first order Ordinary Differential Equation solver. It integrates sets of first order differential equations.
                dy(i)/dt = f(t,y(1),y(2),...,y(N))
                y(i) given at  t .

Syntax: ode ( rhsf, ystart, tstart, tend, dtout, relerr, abserr )

        rhsf    A function that evaluates dy(i)/dt at t. The function
                takes two arguments and returns dy/dt. An example that
                generates dy/dt for Van der Pol's equation is shown
                below. 
Can be user or bltin function.

        ystart  The initial values of y, y(tstart).
Must be column or row vector
        tstart  The initial value of the independent variable.
Can also be a matrix - first element will be used.
        tend    The final value of the independent variable.
Can also be a matrix - first element will be used.
Cannot be same as tstart

        dtout   The output interval. The vector y will be saved at
                tstart, increments of tstart + dtout, and tend. If
                dtout is not specified, then the default is to store
                output at 101 values of the independent variable.

        relerr  The relative error tolerance. Default value is 1.e-6.

        abserr  The absolute error tolerance. At each step, ode
                requires that:

                abs(local error) <= abs(y)*relerr + abserr
                
                For each component of the local error and solution
                vectors. The default value is 1.e-6.


        The Fortran source code for ode() is completely explained and
        documented in the text, "Computer Solution of Ordinary
        Differential Equations: The Initial Value Problem" by 
        L. F. Shampine and  M. K. Gordon.

        Example:
        
        vdpol = function ( t , x ) 
        {
          local (xp)
          xp = zeros(2,1);
          xp[1] = x[1] * (1 - x[2]^2) - x[2];
          xp[2] = x[1];
          return xp;
        };
        
        x0 = [0; 0.25];
        
        xbase = ode( vdpol, x0, 0, 20, 0.05, 1e-9, 1e-9);


Subsections